content
- 線形回帰モデル
- 回帰モデルにおけるベイズ推定
- モデル選択
ベイズ統計学勉強会 夏`22
安藤 正和
⇨所与の年齢と運動療法のもとで酸素摂取量の条件付き分布を推定したい
条件付き分布\(p(y|x)\)が\(x\)の関数として滑らかに変化すると仮定
得られた\(x\)のデータから他の値の情報を得る
条件付き平均\(E[Y|x]\)はパラメータに関して線形であると定める
\[ \int yp(y|\boldsymbol{x})dy = E[Y|\boldsymbol{x}]=\beta_1x_1+...+\beta_px_p=\boldsymbol{\beta}^T\boldsymbol{x} \]
\[ Y_i = \beta_1x_{i,1}+\beta_2x_{i,2}+\beta_3x_{i,3}+\beta_4x_{i,4}+\epsilon_i\tag{9.1} \]
今回のモデルでの\(Y\)の条件付き期待値は、\(x_{i,2}\)のとりうる値によって次のようになる
\[ E[Y|\boldsymbol{x}] = \beta_1+\beta_3\times年齢(x_2=0の場合)\\ E[Y|\boldsymbol{x}] = (\beta_1+\beta_2)+(\beta_3+\beta_4)\times年齢(x_2=1の場合) \]
年齢との線形関係は運動療法のグループ間で切片と傾きの違いがあることが仮定
\[ \epsilon_1,...,\epsilon_n\sim \mathrm{i.i.d. normal}(0,\sigma^2)\\ Y_i=\boldsymbol{\beta}^T\boldsymbol{x}_i+\epsilon_i \]
→\(\boldsymbol{x}_i,\boldsymbol{\beta},\sigma^2\)で条件づけたもとで観測データ\(y_1,...y_n\)の同時分布を完全に特定する
同時確率密度は式(9.2)で書ける
\[ p(y_1,...,y_n|\boldsymbol{x}_1,...,\boldsymbol{x}_n,\boldsymbol{\beta},\sigma^2)\tag{9.2}\\ \]
\[ =\Pi_{i=1}^n p(y_i|\boldsymbol{x}_i,\boldsymbol{\beta},\sigma^2)\\ =(2\pi\sigma^2)^{-n/2}\mathrm{exp}\{-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\boldsymbol{\beta}^T\boldsymbol{x}_i)^2\}\tag{9.3} \]
この同時確率密度は多変量正規分布を用いて書ける
\[ \{\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\sigma^2\}\sim \mathrm{multivariate\ normal}(\boldsymbol{X\beta},\sigma^2\mathrm{\boldsymbol{I}}) \]
\(\boldsymbol{X\beta}\)は以下で示せる
\[ \boldsymbol{X\beta}=\begin{pmatrix}x_1\rightarrow \\ x_2 \rightarrow\\ \vdots \\ x_n \rightarrow\end{pmatrix} \begin{pmatrix}\beta_1 \\ \beta_2\\ \vdots \\ \beta_p\end{pmatrix} =\begin{pmatrix}\beta_1x_{1,1}+\dots\beta_px_{1,p} \\ \vdots \\ \beta_px_{n,1}\dots\beta_px_{n,p}\end{pmatrix} =\begin{pmatrix}\mathrm{E}[Y_1|\boldsymbol{\beta},\boldsymbol{x}_1]\\ \vdots \\ \mathrm{E}[Y_n|\boldsymbol{\beta},\boldsymbol{x}_n]\end{pmatrix} \]
\[ (2\pi\sigma^2)^{-n/2}\mathrm{exp}\{-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\boldsymbol{\beta}^T\boldsymbol{x}_i)^2\}\tag{9.3} \]
式(9.3)の密度は残差\((y_i-\boldsymbol{\beta}^T\boldsymbol{x}_i)\)を通じて\(\boldsymbol{\beta}\)に依存している
観測されたデータを所与とすると、残差平方和\(\mathrm{SSR(\boldsymbol{\beta})}=\sum_{i=1}^n(y_i-\boldsymbol{\beta}^T\boldsymbol{x}_i)^2\)を最小にすることで尤度が最大になる
残差平方和を最小にするには微分する
\[ \mathrm{SSR(\boldsymbol{\beta})}=\sum_{i=1}^n(y_i-\boldsymbol{\beta}^T\boldsymbol{x}_i)^2=(\boldsymbol{y}-\boldsymbol{X\beta})^T(\boldsymbol{y}-\boldsymbol{X\beta})\\ =\boldsymbol{y}^T\boldsymbol{y}-2\boldsymbol{\beta}^T\boldsymbol{X}^T\boldsymbol{y}+{\beta}^T\boldsymbol{X}^T{X}\boldsymbol{\beta} \]
\[ \frac{d}{d\boldsymbol{\beta}}\mathrm{SSR}(\boldsymbol{\beta})=\frac{d}{d\boldsymbol{\beta}}(\boldsymbol{y}^T\boldsymbol{y}-2\boldsymbol{\beta}^T\boldsymbol{X}^T\boldsymbol{y}+{\beta}^T\boldsymbol{X}^T{X}\boldsymbol{\beta})\\ =-2\boldsymbol{X}^T\boldsymbol{y}+-2\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta} \]
\[ \frac{d}{d\boldsymbol{\beta}}\mathrm{SSR}(\boldsymbol{\beta})=0\Leftrightarrow2\boldsymbol{X}^T\boldsymbol{y}+-2\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta}=0\\ \Leftrightarrow \boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta}=\boldsymbol{X}^T\boldsymbol{y}\\ \Leftrightarrow \boldsymbol{\beta}=(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y} \]
式(9.1)のモデルにおいて最小二乗推定量を求め、運動療法の違いを評価する
\[ \boldsymbol{y}= (-0.87,-10.74,-3.27,-1.97,7.50,-7.25,17.05,4.96,10.40,11.05) \]
[,1] [,2] [,3] [,4]
[1,] 1 0 23 0
[2,] 1 0 22 0
[3,] 1 0 22 0
[4,] 1 0 25 0
[5,] 1 0 27 0
[6,] 1 0 20 0
[7,] 1 1 31 31
[8,] 1 1 23 23
[9,] 1 1 27 27
[10,] 1 1 28 28
[11,] 1 1 22 22
[12,] 1 1 24 24
\(\hat{\boldsymbol{\beta}}_{ols}=(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}\)を計算する
[,1]
[1,] -51.2939459
[2,] 13.1070904
[3,] 2.0947027
[4,] -0.3182438
※正しい回帰直線はスライド11の(\(\beta_2\neq0,\ \beta_4\neq0\))
\(\sigma^2\)の不偏推定量は\(\mathrm{SSR}(\hat{\boldsymbol{\beta}}_{ols})/(n-p)\)
\(\boldsymbol{x}_i,\boldsymbol{\beta},\sigma^2\)で条件づけたもとで観測データ\(y_1,...y_n\)の同時分布を式(9.3)に示した
\[ p(y_1,...,y_n|\boldsymbol{x}_1,...,\boldsymbol{x}_n,\boldsymbol{\beta},\sigma^2)\\ =(2\pi\sigma^2)^{-n/2}\mathrm{exp}\{-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\boldsymbol{\beta}^T\boldsymbol{x}_i)^2\}\tag{9.3} \]
このデータ密度は、\(\boldsymbol{\beta}\)の関数として次の様にかける
\[ p(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\sigma^2)\propto \mathrm{exp}\{-\frac{1}{2\sigma^2}\sum_{i=1}^n\mathrm{SSR}(\boldsymbol{\beta})\}\\ =\mathrm{exp}\{-\frac{1}{2\sigma^2}\sum_{i=1}^n [\boldsymbol{y}^T\boldsymbol{y}-2\boldsymbol{\beta}^T\boldsymbol{X}^T\boldsymbol{y}+\boldsymbol{\beta}^T\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta}]\} \]
もし\(\boldsymbol{\beta}\sim \mathrm{multivariate normal}(\boldsymbol{\beta}_0,\boldsymbol{\sum}_0)\)ならば
\[ p(\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^2) \propto p(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\sigma^2)\times p(\boldsymbol{\beta})\\ \propto \mathrm{exp}\{-\frac{1}{2}(-2\boldsymbol{\beta}^T\boldsymbol{X}^T\boldsymbol{y}/\sigma^2+\boldsymbol{\beta}^T\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\beta}/\sigma^2) -\frac{1}{2}(-2\boldsymbol{\beta}^T\sum_0^{-1}\boldsymbol{\beta}+\boldsymbol{\beta}^T\sum_0^{-1}\boldsymbol{\beta})\}\\ =\mathrm{exp}\{ \boldsymbol{\beta}^T(\sum_0^{-1}\boldsymbol{\beta}_0+\boldsymbol{X}^T\boldsymbol{y}/\sigma^2) -\frac{1}{2}\boldsymbol{\beta}^T(\sum_0^{-1}+\boldsymbol{X}^T\boldsymbol{X})/\sigma^2)\boldsymbol{\beta}\} \]
これは、多変量正規分布の密度関数であり、平均と分散は以下で示される
\[ \mathrm{Var}[\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^2]=(\sum_0^{-1}+\boldsymbol{X}^T\boldsymbol{X}/\sigma^2)^{-1}\tag{9.4} \]
\[ \mathrm{E}[\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^2]=(\sum_0^{-1}+\boldsymbol{X}^T\boldsymbol{X}/\sigma^2)^{-1}(\sum_0^{-1}\boldsymbol{\beta}_0+\boldsymbol{X}^T\boldsymbol{y}/\sigma^2)\tag{9.5} \]
ほとんどの正規モデルでは、\(\sigma^2\)の準共役事前分布は逆ガンマ分布である
(5. 正規モデル参照)
\(\gamma=1/\sigma^2\)を観測の精度とすると、\(\gamma\sim\mathrm{gamma}(\nu_0/2,\nu_0\sigma_o^2/2)\)ならば
\[ p(\boldsymbol{\gamma}|\boldsymbol{y},\boldsymbol{X},\boldsymbol{\beta})\propto p(\boldsymbol{\gamma})p(\boldsymbol{y}|\boldsymbol{\gamma},\boldsymbol{X},\boldsymbol{\beta})\\ \propto[\gamma^{\nu_0/2-1}\mathrm{exp}(-\gamma\times\nu_0\sigma_0/2)] \times[\gamma^{n/2}\mathrm{exp}(-\gamma\times\mathrm{SSR}(\boldsymbol{\beta})/2)]\\ =\gamma^{(\nu_0+n)/2-1}\mathrm{exp}(-\gamma[\nu_0\sigma_o^2+\mathrm{SSR}(\boldsymbol{\beta})]/2 \]
これはガンマ分布の密度だとみなせるので、以下の結果が得られる
\[ \{\sigma^2|\boldsymbol{y},\boldsymbol{X},\boldsymbol{\beta}\}\sim \mathrm{inverse\ gamma}([\nu_0+n]/2, [\nu_0\sigma_0^2+\mathrm{SSR}(\boldsymbol{\beta})]/2) \]
※どんな感じで事前分布のパラメータの値を決めてるかはp.173参照
どうする?
弱情報事前分布の例として単位情報事前分布(unit information prior)をあげる
(Kass and Wasserman, 1995)
1個の観測地に含まれる情報と同じ量の情報を含む事前分布
よって、単位情報事前分布では\(\sum_0^{-1}=(\boldsymbol{X}^T\boldsymbol{X})/(n\sigma^2)\)
\(\boldsymbol{\beta}_0=\hat{\boldsymbol{\beta}}_{ols}\)と定めることで、\(\boldsymbol{\beta}\)の事前分布の中心を最小二乗推定量にすることが提案されている
単位情報事前分布は、\(y\)に含まれる情報を少ししか使ってない
同様に\(\nu_o=1,\sigma_o^2=\hat{\sigma}_{ols}^2\)と定めることで、\(\sigma^2\)の事前分布の中心を\(\hat{\sigma}_{ols}^2\)とする
パラメータ推定は説明変数の単位の返還に関して不変であるべき
例: 酸素摂取量を分析するときに用いた年齢(\(x_{i,3}\))が月齢(\(\tilde{x}_{i,3}\))で事後分布は変わらないようにする
与えられた説明変数\(\boldsymbol{X}\)について、ある\(p\times p\)行列\(\boldsymbol{H}\)があって\(\tilde{\boldsymbol{X}}=\boldsymbol{XH}\)
不変原理において、\(y\)と\(\boldsymbol{X}\)から\(\tilde{\boldsymbol{\beta}}\)の事後分布をえたとき、\(\boldsymbol{\beta}=\boldsymbol{H\tilde{\beta}}\)
この条件が成立するのは、任意の正の値\(k\)に対して\(\boldsymbol{\beta}_0=0\)かつ\(\sum_0=k(\boldsymbol{X}^T\boldsymbol{X})^{-1}\)となるとき
\(k\)の特定の仕方として、代表的なものは誤差の分散\(\sigma^2\)に関係させて\(k=g\sigma^2\)とすること
\[ \mathrm{Var}[\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^2] =[\boldsymbol{X}^T\boldsymbol{X}/(g\sigma^2)+\boldsymbol{X}^T\boldsymbol{X}/\sigma^2]^{-1}\\ =\frac{g}{g+1}\sigma^2(\boldsymbol{X}^T\boldsymbol{X})^{-1}\tag{9.6} \]
\[ \mathrm{E}[\boldsymbol{\beta}|\boldsymbol{y},\boldsymbol{X},\sigma^2] =[\boldsymbol{X}^T\boldsymbol{X}/(g\sigma^2)+\boldsymbol{X}^T\boldsymbol{X}/\sigma^2]^{-1}\boldsymbol{X}^T\boldsymbol{y}/\sigma^2\\ =\frac{g}{g+1}\sigma^2(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}\tag{9.7} \]
\[ (\boldsymbol{y}|\boldsymbol{X},\sigma^2) =\int(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\sigma^2)(\boldsymbol{\beta}|\boldsymbol{X},\sigma^2)d\boldsymbol{\beta} \]
積分の中にある二つの密度は式(9.6)を\(V\)と式(9.7)を\(m\)といて以下のように書ける
(式の展開はp.175~176参照)
\[ (\boldsymbol{y}|\boldsymbol{X},\sigma^2) =\int(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\beta},\sigma^2)(\boldsymbol{\beta}|\boldsymbol{X},\sigma^2)d\boldsymbol{\beta}\\ =[(2\pi\sigma^2)^{-n/2}\mathrm{exp}(-\frac{1}{2\sigma^2}\boldsymbol{y}^T\boldsymbol{y})] \times [(1+g)^{-p/2}\mathrm{exp}(\frac{1}{2}m^TV^{-1}m)] \]
ここで、指数の項を組み合わせることで
\[ (\boldsymbol{y}|\boldsymbol{X},\sigma^2)=(2\pi)^{-2/n}(1+g)^{-p/2}(\sigma^2)^{-n/2}\mathrm{exp}(\frac{1}{2}\mathrm{SSR}_g) \]
\[ \mathrm{SSR}_g =\boldsymbol{y}^T\boldsymbol{y}-\sigma^2\boldsymbol{m}^T\boldsymbol{V}^{-1}\boldsymbol{m} =\boldsymbol{y}^T(\mathrm{I}-\frac{g}{g+1}\boldsymbol{X}(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T)\boldsymbol{y} \]
\[ p(\gamma|\boldsymbol{y},\boldsymbol{X})\propto p(\gamma)(\boldsymbol{y}|\boldsymbol{X},\gamma)\\ \propto \mathrm{dgamma}(\gamma,[\nu_0+n]/2,[\nu_0,\sigma_o^2+\mathrm{SSR}_g]/2) \]
\[ \{\sigma^2|\boldsymbol{y},\boldsymbol{X}\}\sim \mathrm{inverse\ gamma}([\nu_0+n]/2,[\nu_0,\sigma_o^2+\mathrm{SSR}_g]/2) \]
まとめると、
intercept aerobic age aerobic.age
-47.592 12.136 1.943 -0.295
#### 事後分布からモンテカルロ標本を生成するコード
yX.o2uptake <-
structure(c(-0.87, -10.74, -3.27, -1.97, 7.5, -7.25, 17.05, 4.96,
10.4, 11.05, 0.26, 2.51, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 23, 22, 22, 25, 27, 20, 31,
23, 27, 28, 22, 24, 0, 0, 0, 0, 0, 0, 31, 23, 27, 28, 22, 24),
.Dim = c(12L, 5L),
.Dimnames = list(NULL, c("uptake", "intercept", "aerobic",
"age", "aerobic.age")))
y <- yX.o2uptake[,1]; X <- yX.o2uptake[,-1] # データ
g <- length(y); nu0 <- 1; s20 <- 8.54 # 事前分布パラメータ
S <- 1000 #サンプリングサイズ
n <- dim(X)[1]; p <- dim(X)[2]
Hg <- (g/(g+1)) * X %*% solve(t(X)%*%X) %*% t(X)
SSRg <- t(y) %*% (diag(1, nrow=n) - Hg) %*% y #残差平方和
s2 <- 1 / rgamma(S, (nu0+n)/2, (nu0*s20+SSRg)/2) #sigma
Vb <- g * solve(t(X)%*%X) / (g+1) #分h散
Eb <- Vb %*% t(X) %*% y #平均
E <- matrix(rnorm(S*p,0,sqrt(s2)), S, p)
beta <- t( t(E%*%chol(Vb))+ c(Eb)) #beta
round( apply(beta,2,mean), 3) Xintercept Xaerobic Xage Xaerobic.age
-47.3482578 12.0988527 1.9335717 -0.2937635
モデルどうする?
説明変数64個
主効果項: 説明変数10個
交互作用項: \(\tbinom{10}{2}=45\)個
二次の項: 9個(性別除く)
変数は標準化(\(y,\boldsymbol{X}\)が平均0,分散1)
[1] 0.669
[1] "平均二乗誤差: 0.668828410609429"
この方法は、変数減少法(backward elimination)の一種である
[1] 0.6594587
仮に、すべての説明変数が\(Y\)と独立な場合、この手順によってどのようなモデルが得られるのか
実際には変数間は独立だが、変数減少法によって誤った示唆がなされてしまうことがある(Berk, 1978)
この時、回帰式は
\[ y_i=z_ib_1x_{i,1}+,\dots,+z_pb_px_{i,p}+\epsilon_i \]
酸素摂取量データでは
\[ E[Y|x,b,z=(1,0,1,0)]=b_1x_1+b_3x_3\\ E[Y|x,b,z=(1,1,0,0)]=b_1x_1+b_2x_2 \]
→g事前分布が使えそう
\[ p(z|y,\boldsymbol{X}) =\frac{p(z)p(y|\boldsymbol{X},z)}{\sum_{\tilde{z}}p(\tilde{z})p(y|\boldsymbol{X},\tilde{z})} \]
任意の二つのモデル比較は、事後オッズを比較することでも可能
\[ \mathrm{odds}(z_a,z_b|y,\boldsymbol{X}) =\frac{p(z_a|y,\boldsymbol{X})}{p(z_b|y,\boldsymbol{X})} =\frac{p(z_a)}{p(z_b)} \times\frac{p(y|\boldsymbol{X},z_a)}{p(y|\boldsymbol{X},z_b)} \]
\[ p(y|\boldsymbol{X},z)=\int\int p(y,\boldsymbol{\beta},\sigma^2|\boldsymbol{X},z)d\beta d\sigma^2\\ =\int\int p(y|\boldsymbol{\beta},\boldsymbol{X},z,\sigma^2)p(\boldsymbol{\beta}|\boldsymbol{X},z,\sigma^2)p(\sigma^2) d\beta d\sigma^2 \tag{9.8} \]
\[ \{\boldsymbol{\beta}_z|\boldsymbol{X}_z,\sigma^2\} \sim\mathrm{multivariate\ normal}(\boldsymbol{0},g\sigma^2[\boldsymbol{X}_z^T\boldsymbol{X}_z]^{-1}) \]
式(9.8)で\(\boldsymbol{\beta}\)を積分
\[ p(y|\boldsymbol{X},z)=\int(\int p(y|\boldsymbol{X},z,\sigma^2,\boldsymbol{\beta})p(\boldsymbol{\beta}|\boldsymbol{X},z,\sigma^2)d\boldsymbol{\beta})p(\sigma^2)d\sigma^2\\ =\int p(y|\boldsymbol{X},z,\sigma^2)p(\sigma^2)d\sigma^2 \]
\[ p(y|\boldsymbol{X},z,\gamma)\times p(\gamma)\\ =(2\pi)^{-n/2}(1+g)^{-p_z/2}\times[\gamma^{n/2}e^{-\gamma\mathrm{SSR}_g^z/2}] \times(\nu_0\sigma_0^2/2)\Gamma(\nu_0/2)^{-1}[\gamma^{\nu_0/2-1}e^{-\gamma\nu_0\sigma_0^2/2}]\tag{9.9} \]
\[ \gamma^{(\nu_0+n)/2-1}\mathrm{exp}(-\gamma\times(\nu_0\sigma^2_0+\mathrm{SSR}_g^z)/2) =\frac{\Gamma([\nu_0+n]/2)}{([\nu_0\sigma^2_0+\mathrm{SSR}_g^z]/2)^{(\nu_0+n)/2-1}} \]
よって\(p(y|\boldsymbol{X},z)\)は
\[ p(y|\boldsymbol{X},z)=\pi^{-n/2}\frac{\Gamma([\nu_0+n]/2)}{\Gamma(\nu_0/2)}(1+g)^{-p_z/2}\frac{(\nu_0\sigma_0^2)^{\nu_0/2}}{(\nu_0\sigma^2_0+\mathrm{SSR}_g^z)^{(\nu_0+n)/2}} \]
任意の二つのモデル\(z_a\)と\(z_b\)の比(ベイズファクター)は
\[ \frac{p(y|\boldsymbol{X},z_a)}{p(y|\boldsymbol{X},z_b)} =(1+n)^{p_{z_b}-p_{z_a}/2}(\frac{s_{z_a}^2}{s_{z_b}^2})^{1/2} \times(\frac{s_{z_b}^2+\mathrm{SSR}_g^{z_b}}{s_{z_a}^2+\mathrm{SSR}_g^{z_a}})^{(n+1)/2} \]
酸素摂取量のデータに対する回帰モデルは以下の通りであった
\[ E[Y_i|\boldsymbol{\beta},x_i]=\beta_1x_{i,1}+\beta_2x_{i,2}+\beta_3x_{i,3}+\beta_4x_{i,4} \]
| \(\boldsymbol{z}\) | モデル | \(log\ p(y|\boldsymbol{X},z)\) | \(p(z|y,\boldsymbol{X})\) |
|---|---|---|---|
| (1,0,0,0) | \(\beta_1\) | -44.33 | 0.00 |
| (1,1,0,0) | \(\beta_1+\beta_2\times グループ_i\) | -42.35 | 0.00 |
| (1,0,1,0) | \(\beta_1+\beta_3\times 年齢_i\) | -37.66 | 0.18 |
| (1,1,1,0) | \(\beta_1+\beta_2\times グループ_i+\beta_3\times 年齢_i\) | -36.42 | 0.63 |
| (1,1,1,1) | \(\beta_1+\beta_2\times グループ_i+\beta_3\times 年齢_i\\+\beta_4\times グループ_i\times 年齢_i\) | -37.60 | 0.19 |
\(\beta\)に対してg事前分布
\(\sigma^2\)に対して単位情報事前分布
5つのモデルの中で最も確率が高い者は\(z=(1,1,1,0)\)
年齢を含む三つのモデルの事後確率を足すとほぼ1になる
運動療法の事後確率は0.82
\[ o_j=\frac{\mathrm{Pr}(z_j=1|y,\boldsymbol{X},z_{-j})}{\mathrm{Pr}(z_j=0|y,\boldsymbol{X},z_{-j})} =\frac{\mathrm{Pr}(z_j=1)}{\mathrm{Pr}(z_j=1)} \times \frac{p(y|\boldsymbol{X},z_{-j},z_j=1)}{p(y|\boldsymbol{X},z_{-j},z_j=0)} \]
## Don't run it again if you've already run it
runmcmc<-!any(dir("source") %in% "diabetesBMA.RData")
if(!runmcmc){ load("source/diabetesBMA.RData") }
if(runmcmc){
BETA<-Z<-matrix(NA,S,p)
z<-rep(1,dim(X)[2] )
lpy.c<-lpy.X(y,X[,z==1,drop=FALSE])
for(s in 1:S)
{
for(j in sample(1:p))
{
zp<-z ; zp[j]<-1-zp[j]
lpy.p<-lpy.X(y,X[,zp==1,drop=FALSE])
r<- (lpy.p - lpy.c)*(-1)^(zp[j]==0)
z[j]<-rbinom(1,1,1/(1+exp(-r)))
if(z[j]==zp[j]) {lpy.c<-lpy.p}
}
beta<-z
if(sum(z)>0){beta[z==1]<-lm.gprior(y,X[,z==1,drop=FALSE],S=1)$beta }
Z[s,]<-z
BETA[s,]<-beta
if(s%%10==0)
{
bpm<-apply(BETA[1:s,],2,mean) ; plot(bpm)
cat(s,mean(z), mean( (y.te-X.te%*%bpm)^2),"\n")
Zcp<- apply(Z[1:s,,drop=FALSE],2,cumsum)/(1:s)
plot(c(1,s),range(Zcp),type="n") ; apply(Zcp,2,lines)
}
}
save(BETA,Z,file="source/diabetesBMA.RData")
} [1] 0.0002349689 -0.0859065610 0.3423203546 0.1784059483 -0.3762771126
[6] 0.2781232114 -0.0044069481 0.0266849669 0.4270552871 0.0056635027
[11] 0.0042001884 0.0146388944 0.0047480761 -0.0021133196 -0.0033337196
[16] -0.0003881057 0.0114732862 0.0132863289 0.0114692165 0.1015032151
[21] 0.0010538499 0.0094550788 -0.0014360247 -0.0241805632 0.0027968474
[26] 0.0025953718 0.0450453591 0.0082762134 0.0149977516 0.0133097496
[31] 0.0015971870 -0.0010502676 0.0025582091 -0.0001147017 0.0010068874
[36] 0.0061732992 0.0159645414 -0.0012871720 0.0003351759 -0.0012023777
[41] 0.0014035074 0.0011455487 0.0201005590 0.0051882624 0.0033907763
[46] 0.0005978504 0.0022609082 0.0111657754 0.0001771299 -0.0005611542
[51] 0.0051321670 -0.0210297438 -0.0059695252 0.0059841046 -0.0059704014
[56] 0.0171164352 0.0053458132 0.0054515426 -0.0031012637 0.0058710738
[61] 0.0008870141 -0.0180022426 0.0044849865 0.0040386379
[1] 0.4825364
ベイズ統計学勉強会 夏`22